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Abstract 

The paper is concerned with dynamics of multi-phase media consisting 
of a solid permeable material and a compressible Newtonian fluid. Gov- 
erning macroscopic equations are derived starting from the space-averaged 
microscopic mass and momentum balances. The Reynolds stress models 
(i.e., momentum dispersive fluxes) are discussed, and a suitable model is 
developed. In the case of granular media the solid constituent is consid- 
ered as an elastic-plastic matrix, and the yield condition is approximated 
by Coulomb friction law. It is revealed that the classical principle of max- 
imum plastic work is not, in general, valid for granular media, and an 
appropriate variational principle is developed. This novel principle coin- 
cides with the maximum plastic work principle for the case of cohesionless 
granular media. 

1 Introduction 

Multi-phase mixtures play a vital part in many natural phenomena and branches 
of engineering (e.g., @], [10], [TS], PU, [22], [25], [27], [25], [25]), and hence the 
development of multi-phase dynamics is of great scientific and industrial impor- 
tance. We restrict our consideration to isotropic (e.g., @], [23], [32]) permeable 
(granular, porous, etc'.) media consisting of a solid matrix and a compressible 
Newtonian fluid. Such media have received the most study (see, e.g., [1], [2], [4], 
[T2] . [27], [28] , [29] , [33], [37]), nevertheless a number of important problems still 
remain to be solved. This paper is mainly concerned with the two important 
problems, namely, modeling of dispersive flux of momentum [1] as well as devel- 
opment of a variational principle for plastic deformations [15] of fluid-saturated 
granular media. 

The so-called macroscopic balance equations are mainly developed by a 
method of averaging (see, e.g., [J], [TT], [27], [23]) of micro-equations. The 
method of averaging over elementary volume of a multi-phase medium contain- 
ing the full ensemble of realizations was originally suggested by Nikolaevskiy 
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et al. (1970) (see the references in [57], [H]). This approach, in contrast to 
the purely phenomenological one (e.g., [16], [28]) gives a possibility to evaluate 
theoretically the type of constitutive laws and sometimes the values of rheologi- 
cal parameters [37]. However, because of non-linearity, the averaging procedure 
leads to the macro-equations that are similar to Reynolds equations (e.g., [3J, 
[23 J for turbulent flows, and hence the problem of Reynolds stress (i.e., disper- 
sive flux of momentum [4]) modeling is coming into play. In the fluid phase the 
pseudo-Brownian velocity pulsations exist even under small Reynolds number as 
a result of chaotic micro-structure of multi-phase media, and these pulsations 
are mainly changing in space [28j . For such a motion the dispersivity tensor 
[35] in the developed models of convective diffusion in, e.g., porous media [57J 
sec. 7.3] depends on the average velocity, whereas, because of Galilei-Newton 
principle of relativity [38] V. 2, p. 105], the turbulence model parameters do 
not depend on the velocity as such (see, e.g., [23], [3]). Hence, the models of 
turbulent diffusion, as they exist, are not, in general, applicable to multi-phase 
media. In the case of granular media the Reynolds stresses (i.e., momentum 
dispersive fluxes) are often believed as negligible 27] . However, in develop- 
ing models for correct simulation of processes accompanying, e.g., fluid flow in 
vicinities of gas wells having high production, underground explosion works, oil 
and water wells rehabilitation and stimulation by shock technologies (steam in- 
jection, aggressive pressure pulsing, etc.'), meteorological flows over urban and 
vegetative canopies, and so on, the momentum dispersive fluxes must be taken 
into consideration (e.g., [4], [10], [18], [28] V 

A dispersive flux model of an extensive quantity was suggested in [4] for 
a microscopically laminar flow regime. The model is essentially based on the 
modified rule [U Eq. 2.3.48] for volume averaging of a spatial derivative. Since 
the modified averaging rule is also used in the development of the total viscous 
resistance expression as well as the macroscopic momentum balance equation 
for a fluid phase [4] sec. 2.6], we consider the applicability of this rule to a 
viscous fluid flow through a permeable medium. 

The authors [I] attempted to develop the modified rule for a quantity, G, that 
attains no maximum or minimum value within the void space of a representative 
elementary volume (REV). On this basis it is assumed [U p. 125] that the 
quantity G is a harmonic function on the microscopic level. This demand is 
sufficient but not necessary condition for the validity of the maximum principle, 
and hence it could be too restrictive. To demonstrate it we assume, for the 
sake of simplicity, that the solid matrix is immobile and the entire void space is 
occupied by a single Newtonian fluid with the density p — const and the dynamic 
viscosity u = const. Eliminating the body force of gravity by subtraction from 
the true pressure p of the hydrostatic pressure, the Navier-Stokes and continuity 
equations (e.g., [3J, [23] . [33]) can be written in the following non-dimensional 
form: 

5 ^ + ( V ' V)V = -£ U V?+^V 2 V, (1) 
V-V = 0, S h =-^-,E u = ^,Re=^, (2) 
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where V, P are the non-dimensional velocity and pressure, respectively; v de- 
notes the kinematic viscosity; Sh, E u , Re are, respectively, Struhal, Euler, and 
Reynolds numbers; the reference quantities are denoted by an asterisk. Multi- 
plying ((T|) by V we obtain, in view of the first equation in that 



Thus, the basic equality, i.e., the Laplace equation for the pressure [4] Eq. 
2.6.6], can be approximately valid if the right-hand side in ((3]) is negligible, i.e., 
in general, if E u ^> 1. Hence, we conclude that the derivation of the viscous 
resistance expression as well as the macroscopic momentum balance equation 
[U sec. 2.6], relying on the modified averaging rule, can be, in general, valid for 
a creeping flow only. 

The applicability of the modified averaging rule for the development of the 
momentum dispersive flux model is based on the assumption that the linear 
momentum density (pV) is a harmonic function on the microscopic level. In 
view of ([1]) we obtain that pV will be an approximate solution of the Laplace 
equation if the following inequalities will be valid simultaneously: ShRe <C 1, 
Re <C 1, E u Re « 1. In general, it is possible for a very specific flow. As an 
example, let us consider the Newtonian fluid flow (see the governing equations in 
(fTJ) and |(5J)) through a porous medium as made up of a bundle of parallel tubes 
whose radii are assumed to be uniform in size. In the case of a steady-state 
laminar flow regime we have (in every tube), in fact, Hagen-Poiseuille flow [351 
p. 117]. In polar coordinates the velocity (V) distribution over a cross-section 
of a tube can be written in the form 



where R denotes the radius of the tube; L denotes the full length of the tube; 
C denotes the pressure drop; Po, Pl are the pressures at the bases of the tube. 
Obviously, the linear momentum density (pV) in the case of Hagen-Poiseuille 
flow is not a harmonic function. To elucidate the conditions wherein pV is 
an approximate solution of the Laplace equation we assume that the reference 
velocity is equal to the maximum value of the velocity in Hagen-Poiseuille flow 
(i.e., V* — 0.25C-R 2 / p), I* — R, and P* = Po — Pl- In such a case we obtain 
that E u Re = AL/R, and hence pV will be an approximate solution of the 
Laplace equation if R 3> L. The last inequality is to say that the influence 
of the boundary (non-slip) conditions on the flow regime must be negligible. 
This is by no means the case of porous media. Thus, we conclude that the 
applicability of the modified averaging rule to the linear momentum density 
(pV) is, in general, questionable. 

Using the modified averaging rule as well as a number of additional assump- 
tions and approximations, Bear and Bachmat [4] found that the dispersive flux 
of an extensive quantity is proportional to the gradient of the mean density of 
the extensive quantity. Hence, the momentum dispersive flux is proportional to 



V 2 P 



— V-[(V-V)V]. 



(3) 
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(4) 
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the gradient of the mean momentum density of mass. The mean momentum 
density in its turn may be decomposed into two fluxes: a macroscopic advective 
flux and a dispersive flux of mass. The latter flux, in view of the model by 
Bear and Bachmat [3], is proportional to the gradient of the mean density of 
mass. Thus, following [3], we obtain an interesting result that the macroscopic 
momentum balance equation is a third-order partial differential equation. How- 
ever, in view of the foregoing analysis, the basis for this result is not convincing. 
Two more points need to be made. The momentum dispersive flux (see, e.g., 
left-hand side in [JJ Eq. 2.6.52]), which is a symmetrical tensor, is approximated 
by the tensor (see the right-hand side in [H Eq. 2.6.52]), which is not, in gen- 
eral, symmetric. Furthermore, the kinetic energy of dispersion (analog to the 
kinetic energy of turbulence, see [3j p. 221]) is not taken into consideration in 
the momentum dispersive flux model in [4]. 

Let us note, however, that the approach employed by Bear and Bachmat [4] 
deserves more attention. This approach and the conclusions, as it is underlined 
by the authors, are equally valid for any extensive quantity. It is a reflection of 
the plausible assumption that the mechanisms of heat, mass, and momentum 
transfer in multi-phase media are identical. In the semi-empirical theory of tur- 
bulence the similar assumption is called as Reynolds analogy (see, e.g., [3J, [2"31 . 
[55]). Specifically, Anderson et al. [3J pointed out that the ratio of the diffu- 
sivities for the turbulent transport of heat and momentum (turbulent Prandtl 
number) is a well-behaved function across the flow, and the Prandtl number 
varies between about 0.6 at the outer edge of the boundary layer to about 1.5 
near the wall. Assuming that this assumption is valid, we can conclude [4] that 
the coefficient of mechanical dispersion (dispersivity tensor, see [32]) has anal- 
ogous form for any conservative extensive quantity, be it heat, mass, or linear 
momentum. 

It should also be remarked that Bear and Bachmat [4] suggested the totally 
irreversible model of the momentum dispersive flux for the case of microscop- 
ically laminar flow regime provided the density is constant. In contrast to [3], 
a reversible model of the Reynolds stress tensor (i.e., the momentum dispersive 
flux) is developed by Nikolaevskiy [2S]. It is assumed that the local velocity 
(V) is a stochastic function of the average velocity (V). In such a case the 
vector of mean velocity is transformed randomly (see [26], [28]) into the velocity 
pulsations (V = V— V). This stochastic transformation is represented, in both 
the symbolic and indicial notations, as follows: 

V' = L V; Vl = UjVi, (5) 

where the repeated indices, as usual, denote summation; the tensor L is de- 
termined by the structure of porous medium, Re, and a realization parameter, 
corresponding to the random character of the medium. Then, the equations for 
the Reynolds stresses are, in fact, written 28, sec. 4.4.1] in the following form 

V¥ = T : (V V) ; VjVj = T ijkl V k V l , T ijkl = , (6) 

To estimate the kinetic energy of dispersion, 0.5K'V^', we assume, in view of 
the isotropy of the porous medium [26j , that the second rank tensor with the 
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components L^Lu will be isotropic. Then, by virtue of ([S]), we obtain 



{\'f = VlV( = trVlV; = L ik LaV k Vi=w(yf, u = -L : L. (7) 

Let us assess the validity of the assumption fl5j and, hence, the model ([6]). 
Notice, if the direction of V in ([5j will be changed to the opposite one (i.e., V — > 
— V), then the only direction of V will be changed (i.e., V — > —V). However, 
such a property is not, in general, exhibited by the flow, e.g., in convergent 
and divergent channels. The exact solution of the Navier-Stokes equations for 
such a flow was originally found by Jeffery and Hamel (see the brief sketch in 
[35l pp. 104-106]). The velocity distribution for the convergent and for the 
divergent channel differ significantly from each other, and in the latter case 
vary greatly with Reynolds number (see, e.g., |13j . |35|). Inasmuch as the flow 
in such channels is essentially irreversible, we conclude that ^ can, in general, 
be approximately valid under low Reynolds numbers only. 

Frick [13) points out that a symmetric divergent Jeffery-Hamel flow exists 
only if the Reynolds number Re < Re and the opening angle a < a, where the 
values Re and a meet the following condition: Re — 6 (ir 2 /a — a) . Therefore, 
choosing a > tt we obtain at least one region of back-flow, whichever Re might 
be, and, resulting from it, separation. Separation of a boundary layer, in reality, 
gives rise to vortices [35J Sec. 2] resulting in turbulence. For instance, measure- 
ments demonstrate (351 Sec. 7.2.6] that the flow of a free jet can be laminar 
up until about Re = 30, where Re is referred to the outlet velocity and the slit 
height. Hence, in a real granular medium, the solid phase of which composed, in 
general, of irregular in size and shape grains, the vortices (i.e., micro- vortices) 
can arise under low Reynolds numbers as the result of steep rise in pressure 
at sharp edges, fractures, sudden expansions, etc'. The higher Re, the more 
micro-vortices arise within an REV. Nevertheless, the flow is still laminar. Let 
us note, however, that the vortices are the main source of turbulence [6]. The 
turbulent flow stems from the lose of vortex stability and degradation of vor- 
tex structure on further rise in Re. If Re is over a critical Reynolds number, 
then the vortex structure breaks into turbulence within a part of the REV. The 
higher Re, the most part of the flow will be turbulent. The above speculation is 
supported by experimental results (see, e.g., [5], [7], [21] and references therein), 
and hence it might be differentiated the following regimes of flow: 1) Laminar 
regime, where the resistance to the flow is directly proportional to the mean ve- 
locity (Darcy linear law). 2) Laminar regime, where the resistance is nonlinear 
(Darcy-Hazen-Dupuit-Forchheimer law). 3) Transition regime. 4) Turbulent 
flow. Thus, from the preceding, it appears that the assumption ((5J and, hence, 
the model ^ can be approximately valid in the case of the first flow regime 
(laminar regime, Darcy linear flow), otherwise the validity of (J5j is questionable. 
The same conclusion is valid for the model developed by Bear and Bachmat [3] . 

Recently, several turbulence models have been established for turbulent flows 
(i.e., for the fourth regime) in permeable (granular, porous, etc'.) media (see, 
e.g., [10], [18], [22], [24], [31], [39], [40], and references therein). Let us note that 
the widely used one- and two-equation turbulence models (e.g., [10], [18], [22], 
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EH, EH) [10]) are only valid in the fully turbulent regime, i.e., these models 
are not appropriate for the near- wall region [3] pp. 231-233]. Owing to this, 
the validity of these models for turbulent flows in porous media is questionable. 
Thus, as already noted in [HO p. 63], modeling of the Reynolds stresses for 
multi-phase flows is still at its infancy. 

The present study (Sec. [5J is devoted to the development of a sufficiently 
simple mathematical model for simulating flows in permeable media under all 
regimes. To avoid the specific problems of dispersive flux modeling associated 
with inconstant phase densities (e.g., [3J, U), we will develop a Reynolds form 
of the balance equations in mass-averaged variables [3J p. 201]. Then we will 
develop a model of the dispersive flux assuming that: (i) the Reynolds stresses 
are linearly dependent on the mean velocity derivatives (see, e.g., [3], [4], [19] . 
[23] . [28]), (ii) the Reynolds stress tensor will be isotropic if the mean velocity 
is constant [19], (iii) the Reynolds analogy is valid for multi-phase media. 

Considering the granular medium as an elastic-plastic one (e.g., [16], [T7] . 
[55] , [57] , [55J ) , we are concerned with validity of the classical principle of max- 
imum plastic work, as applied to granular media. It can be seen from (e.g., 
[T5I p. 58], [38j V. 4, ch. 12]) that the associative flow rule (the normality law 
[151 p. 58]) is the necessary condition for validity of the maximum plastic work 
principle. Nikolaevskiy [57], [55] has revealed that in the case of granular me- 
dia the irreversible strain-rate must be determined by the non-associative flow 
rule. Thus, we may conclude that the principle of maximum plastic work |15] is 
not in general valid for granular media. Importance of variational principles in 
physics, including the maximum plastic work principle developed by von Mises, 
Taylor, and Bishop and Hill, is well known (e.g., [14], [15], [34], [34]). In par- 
ticular, Han and Reddy [15] wrote, "The principle of maximum plastic work is 
a vital constituent of the theory of plasticity." Moreover, currently, numerical 
methods for solving the problems of elastic-plastic deformations, including con- 
struction of monotone (e.g., [8], [9]) difference schemes, are based on non-classic 
formulation of variational principles, namely variational inequalities (|15j. [30] . 
[34] ) . Hence, there is a need to develop a proper variational principle for plastic 
deformation of granular media. With this in mind we will assume the validity of 
the Coulomb yield condition (e.g., [55], [57], [55]) and the non-associative flow 
rule (e.g., [27], [28]) . Then the desired variational principle will be rigorously 
deduced on the basis of irreversible thermodynamics (e.g., [14], [27]), as applied 
to granular media. Let us note that in such a development one would like to 
use a theoretical premise instead of the Coulomb friction law that is nothing 
more than an empirical relationship [2 5) . Currently another approach, free of 
Coulomb condition, is suggested by Jiang and Liu [17] . The novel elastic the- 
ory [17j accounts for mechanical yield by a feature of non-linear elasticity only. 
However, from physical point of view, such an approach is not well founded, as 
it is not obvious that solid friction at inter-particle contacts can be totally ac- 
counted for by non-linear elasticity. Furthermore, Jiang and Liu [17j developed 
their theory for the case of cohesionless granular media only. Hence, the use of 
this theory [17], as it exists, for development of the variational principle is out 
of question. 
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2 Momentum transport 



We start our investigation with the volume averaged balance equations [4] . In 
the absence of phase transition, the macroscopic mass and, respectively, mo- 
mentum balance equations for a fluid phase can be written in the form 

|(#7) = -v-(^7), (8) 

d_ 

Of 



Pf y f) =-V-^(p / V / V/-ff7j 



+ Uf J <Tf ' nds + < ^/ g/ ' ( 9 ) 

Sfs 

where pf denotes the fluid density, <j> denotes the void fraction, Vf denotes 
the velocity of the fluid, cr j (= — pfl + t/) denotes the fluid stress tensor, 
Pf denotes the pressure, I denotes the identity tensor, t/ denotes the viscous 
stress tensor, gf denotes the body force vector, U / denotes the volume occupied 
by the fluid phase within an REV, Sf a denotes the surface of the fluid-solid 
interface, n denotes the outward unit vector to Uf on Sf s . Let U s denote the 
volume occupied by the solid phase within the REV, and let e/ and e s denote 
any variables referred to the fluid and solid phases, respectively. Hereinafter ej 
and ej denote the volume averages of e/ and e s over, respectively, Uf and U s . 
Let e} (= pf ef/~p~j) and e a (= p s e s /~p^) denote the mass-averaged variables. 
Following Nikolaevskiy [27] we write the phase interaction force (F) in the form 

F = — / cr / -nds-pV0 = -R-U, U = V/-V^, (10) 
Uf J 

Sfs 

where R (=p,tp(f>'K~ 1 ) denotes the symmetric resistance tensor [27], K denotes 

the tensor of absolute permeability, ip — ijj(Re), Re (= U l*/v) denotes the 

local Reynolds number, v — P-/ Pf denotes the kinematic viscosity, Z* denotes 
the length parameter characterizing the void space. 

Using the mass-averaged variables and by virtue of (fiT)]) we obtain from ([8]) , 
J9J the following balance equations: 

|(0p7) = -v-(0pjv / ) ! (ii) 

+ 0pjgf + V • cjrrj - V • <t>J]V"fV"f, (12) 

where e'j = e/ — denotes the fluctuations of a variable e/. We will also use 
e' f = ef — ej. Since the fluid phase is a Newtonian liquid, we can write (e.g., 
[38] ) that 

r / = 2 M S- K V-V / I, S = i [W/ + (W/)*] , (13) 



- V/ = -<^VP7-R- (V/-V 
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where k (= — () denotes the coefficient of bulk viscosity, ( denotes the second 
coefficient of viscosity, ( )* denotes a conjugate tensor. Estimating the mean 
value of the viscous tensor (TJ) we note j3] that in practice the viscous terms 
involving = —p'jW'^/'pj can be neglected, and hence we obtain 



tJk, 2/JS-kV- V/I, S = -VV / +(VV 



(14) 



Notice, the first equality in (fl~4"l) may be considered as strict equation introducing 
new characteristics ~p and k instead of the mean value of conventional coefficients 
of viscosity. In such a case we obtain that v — ~p,/~pj. 

The dispersion in the case of isotropic granular media is, in general, non- 
isotropic [4] , which is to say that at every point of the flow there must be defined 
a symmetric tensor (D/) of second rank (e.g., [4], [23], [27], [32]) such that the 
dispersive flux of the momentum will be dependent on D/ ■ VVj. Since the 
Reynolds stresses form a symmetric tensor, it is natural to consider this tensor 
as a linear function of the following symmetric one 



D/ • VV 



/ 



D t ■ V V . 



(15) 



We will assume that because of chaotic micro-structure of the isotropic granular 
medium, the Reynolds stress tensor will be isotropic if 3?/ =0. Thus, we may 
write that 

-pjvfv'} = pj(* f + dl), (16) 

where d is a scalar parameter. In view of the last equality in (|10p we obtain 
that Vy = U". Then, assuming that 



|U| 2 «(1+ W/ ) 



u 



(17) 



where ujf (> 0) is a new non-dimensional parameter, we can write 



(v^') = UTJ-2U-U+U • U = (U) 2 - 



U 



U/ U 



{ v < V-) . (18) 

It is often assumed (see, e.g., g], [37]) that (U) 2 » flj] (or (U) 2 w (U) 2 ), and 

hence ujf ~ 0. Notice, such an assumption is not valid for laminar regimes, how- 
ever, uj f can be close to zero for the case of turbulent regime under sufficiently 
high Reynolds number (see Sec. [5]). 

Equating linear invariants of the tensors in (|16|) we obtain, in view of (fl8|) , 
the following estimation of the Reynolds stress tensor 



(19) 
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n, 



3 I-/ 



Vf-V, 



(20) 



Using the Reynolds analogy (e.g., [3], [13]), as applied to multi-phase media, 
and following [55] and [35] we find, for the case of isotropic granular media, that 
the dispersivity tensor Dj can be approximated as follows: 



D 



V ~ ~ 

= F 1 vS ij +F 2 ^U l U ] 



(21) 



where 



Fi = 



a° f Re 2 
1 + ctfRe' 



F 2 = 



1 + (3 f Re 



Re = 



U 



(22) 



a°p bf, Oif, and j3f denote the dimensionless parameters describing the geometry 
of the void space, U l = Vf ;i - V s;i . We find, by virtue of (J2TJ)- (J22J) , that 



D 



f;ij 



U 



u 



Re 



1 + P f Re 



(23) 



where a/ = a* (l + /3ji?e) / (1 + otfRe). The similar approximation of the 



dispersivity tensor is suggested in j4, p. 218]. Notice, the approximation (|23[> 
of the dispersivity tensor D f does not contradict with Galilei-Newton principle 
of relativity (38] V. 2, p. 105], since D/ depends on U = V/ — V s . 

Similarly we obtain the macroscopic balance equations for the solid phase: 



(1 



V s = V-«f'-(l-0)Vp7+(l-#75g.+R-(V/-V.) - ( 25 ) 



= -V • 



(24) 



rs dt 

where a? = (1 — 0) (7f s +pjT) is the Terzaghi effective stress [27]. Notice, fol- 
lowing Nikolaevskiy [57], the Reynolds stress tensor (i.e., momentum dispersive 

flux [?]), (1 — 0)p7V"V", is assumed as negligible for the solid phase. In such 
a case we obtain the conventional mathematical model for elastic-plastic defor- 
mations of granular media, i.e., the system of hyperbolic equations. However, 
as it can be concluded, e.g., from [4], the dispersion of the momentum should, 
in general, be taken into consideration if irreversible deformations take place. 
Then, in view of the Reynolds analogy, the momentum dispersive flux for the 
solid phase can be estimated by application of (fT§j). (f2TJ)) . (fl"5"]) . and ([53]) with 
obvious modifications. In such a case we obtain the system of partial differential 
equations of parabolic type as a mathematical model of multi-phase dynamics. 



3 Maximum principle 

Hereinafter, the considered quantities, such as stress, density, etc.', will be of 
the average ones only. Hence, the sign to indicate the fact of averaging will be 
deleted. 
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To develop a variational principle for plastic deformation of granular media 
we start with Coulomb friction law, as applied to a two-phase granular medium. 
According to Terzaghi's principle [37] the yield condition is formulated to the 
effective stresses: 

C a = -|= lovl + aa f - Y = 0, a f = hrcr f , (26) 
V3 3 

where a T denotes the shear stress intensity in the solid matrix, a (= a (x) > 0) 
denotes the internal friction coefficient, x denotes a hardening parameter, and 
Y = Y (x) denotes the cohesion. Strain increment de of the matrix can be 
divided (e.g., [15], [27]) into elastic (de e ) and plastic (de p ) parts: de = de e +de' p . 
The plastic strains are determined by non-associative flow rule [37] that can be 
written in the form 



tr f + - AVI— | 1 + -Aa ) a f I 



3 V 3 



A, (27) 



where A = A (x) is the dilatancy rate, e p = de p /dt, A = d\/dt. The scalar 
function A = if C a < 0. Following Sedov [Ml V. 4, pp. 145-147], in view of 
the first and second law of thermodynamics, as applied to the solid matrix of 
granular media, we obtain 

(1 ^ )fcT f = ^ + 9:eP -°' (28) 

where 9 — cr^ — (1 — 0) p s dF/de p , S denotes the specific entropy, diS denotes 
the specific entropy variation due to irreversible processes, T denotes the abso- 
lute temperature of the solid phase, F = U — TS = F {e e 1 e p ,T) denotes the 
specific free energy, U denotes the internal energy, q is the heat flux. Thus, the 
energy dissipation is determined by thermodynamic currents q, e p and conju- 
gated forces, according to the bilinear form of (|28[) . In view of Curie principle 
[Tl] , the value of energy dissipation associated with plastic deformation D p = 
9 : e p > 0. Hence, the non-associative flow rule (|2"7|) is bound to be equivalent 
to the following constitutive equation [T4 : e p = L : 6, where L is a fourth-rank 
tensor. It is possible on condition that 

Thus, in the case of granular media, the specific free energy F is a function of 
the first invariant of plastic strain tensor e p . 

Following Sadovskii [53] Section 1.2], we assume that there exists a con- 
vex, possibly non-differentiable, function B (e p ) such that 8 € dB(e p ), where 
dB (e p ) denotes the sub-differential ([30]) [33]) of the function B. It is as- 
sumed [331 Section 1.2] that B (e p ) is a positive-homogeneous function of e p , 
i.e., B (be p ) = bB (e p ) for b > 0. The above sub-differential relationship is 
equivalent ([30], [33]) to 

9:e p -B = max (9 : e p - B*) . (30) 
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The right-hand side of ([3H)l is the Young's transformation [31] of B, and, in view 
of the positive-homogeneousness, is the characteristic function of the convex set 
T = {8 | C a < 0}, i.e. the right-hand side of (j3"0]) is equal to zero for 8 e Y, 
and hence B — 8 : e p , i.e. B = D p . Using the Young's transformation for the 
characteristic function of the convex set T, we obtain 

D p = max ( 8« : e p ) . (31) 
e,er 

In view of pip we obtain the desired variational principle: 

(8-8*):de p >0, 8,8, € T, (32) 

where 8 and 8, denote the tensors associated with actual and, respectively, 
virtual energy dissipation. 



4 Concluding remarcs 

To fulfill the basic laws of conservation by the mean quantities, the averages were 
introduced by different methods. In particular, the velocities and body forces 
are mass-averaged, whereas the densities and stresses are treated as volume 
averaged. Owing to this approach we succeeded in modeling of the Reynolds 
stress tensors under variable densities. Thus, the developed macroscopic balance 
equations are applicable for modeling flows through a permeable medium under 
a wide range of velocities. 

The novel variational principle (|32p , developed for plastic deformation of 
granular media, can be interpreted as the maximum plastic energy dissipation 
principle. If a granular medium is idealized to be cohesionlcss (Y = 0), then in 
view of (f2"5|) 8 = a? , 8, — <x{, and hence the novel variational principle (|32[) 
coincides with the maximum plastic work principle |15j . 



5 Appendix 

Let us estimate the value of uj f in (|17[) by considering the motion of an incom- 
pressible viscous fluid through a medium consisting of straight cylindrical tubes. 
Let us note that because of incompressibility we have e/ = ej , V e/ . 

First we consider Hagen-Poiseuille flow (see, e.g., [35j P- 117]), i.e., a steady- 
state laminar flow regime. In polar coordinates the velocity distribution over a 
cross-section of i — th tube can be written in the form 

C-R 2 ( r 2 \ 

where Ri denotes the radius of i — th tube, Ci denotes the pressure drop, and 
k; = (k x ^ , k y _i, k z ^) denotes the unit vector parallel to the axis of the tube i. 
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By virtue of ([55)1 we obtain that 

R , 



Then clearly 

EJ» ECvg| 

Let ^ = Ci-Rf- Then, in view of (f3"5)) . we find 



U = ^2 = tt¥&*<- (35) 



Cu) 2 = E £»fc^) 2 + (E C»fcy,Q 2 + (E £^,*) 2 / 36 n 
1 ) 6V( EjR f) 2 

By virtue of the elementary inequality, a 2 + b 2 > 2ab, it is an easy matter to 
prove that the numerator in the right-hand side of (j3"6")l is bounded above by 
(ECi) 2 - Hence we get 

(U) 2 < & C * R ^ (37) 
In a similar manner, we obtain 

2 _ , k^R! 



Then obviously 



J* = 2tt I (UOVidr, = -g-gL. (38) 



E^f 48M a £*T 



(U) 2 = ^52 = (39) 



In view of and (j3"T|) we obtain from (|17p that 

w > 4EflfEcy _ L (40) 
3(EC^f) 

The lower bound of uj f follows from Cauchy-Schwarz inequality 

(E^) 2 <E a2 E^- ( 41 ) 



Setting di = Ri and = Ci^li we obtain from (|40p . in view of (j4Tj) . that 
> 1/3- Notice, lo f = 1/3 if the tubes are identical in diameter and their axes 
are parallel to each other. 

Next we consider a turbulent regime. The velocity distribution over a cross- 
section of i-th tube can be written (e.g., [361 p. 563]) in the following form 



U, = U mA ( 1 - ) k 2 < n < Ri, n = const, (42) 
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where U m ^ denotes the maximum velocity. Notice, we assume n = const in 
the empirical formula (|42|) . as the exponent n varies slightly with the Reynolds 
number Re [351 p. 563]. In the perfect analogy to the deduction of (|3"T)). (|3"9"1) , 
and pi)]) , we obtain in the case of turbulent regime that 

1 j " (n+l) 2 (2n+l)2(E^ 2 ) 2 ' (n+l)(n + 2)£i2?' 1 j 

u f >- — n ^ t; — 9 m,t - - 1. 44 

4n2(n + 2)(Ef/ m , l i?. 2 ) 

Let a,i = Ri and hi — U m ,iR\ in (|4"Tj) . Then, in view of (|44l . we obtain 

W ' * 4n>(n + 2)- (45) 

In particular, if n = 6, i.e. i?e = 4 • 10 3 36, p. 563], then we find, by virtue of 
(|4"5"1) , that u)f > 0.027. Thus, in the case of turbulent regime the minimal value 
of ujf is far less than in the case of laminar regime. 
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